Lu Xiaoyong, Zhang Xiaozhang, Zhang Zhizhong. Numerical simulation of metal evaporation based on the kinetic model equation and the direct simulation Monte Carlo method. Chinese Physics B, 2018, 27(12): 124702
Permissions
Numerical simulation of metal evaporation based on the kinetic model equation and the direct simulation Monte Carlo method
Lu Xiaoyong1, 2, Zhang Xiaozhang1, †, Zhang Zhizhong2
Department of Engineering Physics, Tsinghua University, Beijing 100084, China
Research Institute of Physics and Chemical Engineering of Nuclear Industry, Tianjin 300180, China
Metal evaporation on the basis of the kinetic model equations (BGK and S-model) and the direct simulation Monte Carlo (DSMC) method was investigated computationally under the circumstances of collimators existing or not. Numerical data of distributions of number density, bulk velocity and temperature were reported over a wide range of evaporation rate. It was shown that these results reached a good agreement for the case of small evaporation rate, while the deviations became increasingly obvious with the increase of evaporation rate, especially when the collimators existed. Moreover, the deposition thickness over substrate obtained from the kinetic model equations were inaccurate even though the evaporation rate was small. All of the comparisons showed the reliability of the kinetic model equations, which require less computational cost at small evaporation rate and simple structure.
The flow of metal vapor is encountered in many project fields and science research, such as in supersonic molecular beam technology,[1] enrichment of uranium,[2–4] electron beam physical vapor deposition,[5–7] and vacuum technology.[8–10] Measurement of the flow field and macroscopic parameters for metal vapor is very difficult because of many uncertain sources of error and equipment’s thermal stability at high temperature. As a result, obtaining and optimizing the flow characteristics theoretically have attracted many scientists’ attention. For example, Karol Waichman used the iterative method to research the effects of boundary geometry on rarefied condensable vapor flow.[1] In his work, the nonstationary Bhatnagar–Gross–Krool (BGK) equations were calculated for a steady-state solution. Nevertheless, the influence of temperature on collision frequency was ignored. Fan et al. performed calculations of electron beam physical vapor deposition using the direct simulation Monte Carlo (DSMC) method[5] and the results were compared with those obtained from an experiment in which the atomic absorption spectra and deposited film thickness were measured. Calculations and experiments were in perfect agreement with each other at two different evaporation rates. In addition, there are some other kinetic methods, namely the S-model[11] and the McCormack model method,[12] which can also be used to describe this problem. However, the Prandtl number of monatomic gas is not the actual number 2/3 in the BGK model equation. The S-model equation has not been proved to obey the H-theorem until now, although it can provide the correct Prandtl number. Therefore, the applicability of the kinetic model equations is open to question in specific conditions.
Felix Sharipov and his team have developed the confirmation works of the kinetic model equations in rarefied gas flow through a slit at arbitrary pressure ratio. The flow properties were compared with those obtained from the linearized Boltzmann equation and the DSMC method.[8,10,13–16] The results showed that the kinetic model equations could provide an accurate flow field of rarefied gas with a less computational effort. While in the case of metal evaporation, the velocity distribution of metal vapor is strongly nonequilibrium in an evaporation chamber and the boundary conditions are quite different. Although the results of kinetic model equations and the DSMC method correspond with each other qualitatively from previous works, scientists have not carried out a research about the comparison between kinetic model equation and the Monte Carlo method on metal evaporation quantitatively.
The aim of this paper is to calculate the flow field and macroscopic parameters of metal evaporation at different evaporation rates and evaporation chamber structures based on the kinetic model equations, including the BGK model equation and the S-model equation. Then the results are compared with those obtained from the DSMC method, which describes rarefied gas flow successfully, to indicate the deviation between these methods quantitatively. The analysis of the deviation and the discussion of the kinetic model equations’ applicability provide a basis for understanding and ameliorating the experiment process. Moreover, the results can be regarded as references to other analogous rarefied gas flows.
2. Statement of the problem
A two-dimensional (2D) evaporation system, which consists of evaporating source and collimators, is enclosed by several walls as shown in Fig. 1. All of the dimensions of the evaporation chamber and collimators are marked in the structure and the value of them are presented in Table 1. A high-energy electron beam heats the metal ingot, creating a liquid metal pool from which surface metal atoms vaporize. The vapor diffuses rapidly and deposits on the chamber boundary finally. The vapor leaves the source surface with the saturated vapor pressure Ps (Ts), where Ts is the temperature of pool surface. Therefore, the number density of metal vapor at the source surface can be obtained nS = Ps (Ts)/k · Ts, where k is the Boltzmann constant. The velocity obeys the Maxwell velocity distribution f(v) corresponding to the characteristic temperature Ts while only that upward is logical. Each wall in the evaporation chamber has its own temperature TW and absorption rate αW. That is to say, when metal atoms collide with the surface of the wall, a fraction (1 − αW) of them reflect diffusely with the temperature TW and the rest condense on the surface. Because the whole structure of evaporation chamber is z-axis symmetric, calculations are carried out for the right half chamber only and the flow properties of the left half can be obtained by mirroring about the z axis.
Fig. 1. (color online) Schematic diagram of two-dimensional metal evaporation.
In this article, two chamber structures are calculated to confirm the applicability of the kinetic model equations. One is the ideal condition that the collimators are not put into the chamber and all of the walls in the chamber are assumed to be sticking enough to absorb the metal atoms completely. The other is that the collimators are in consideration. The surfaces below the dashed line in Fig. 1 have the same temperatures TW and absorption rates αW which are less than 1, while the absorption rate of the surfaces above the dashed line is equal to 1. In addition, metal atoms that return back to the liquid metal pool dissolve entirely in both two conditions. Obviously, the second one is more accordant with the practical circumstances.
Table 1.
Table 1.
Table 1.
Dimensions of the evaporation domain.
.
Length
H
h
L
l
D
a
Value/cm
30
20
32
4
2
1
Table 1.
Dimensions of the evaporation domain.
.
3. Methods
The 2D metal vapor flow without external force can be described by the stationary kinetic model equation as follows:
where f(x, z, vx, vy, vz) is the velocity distribution function and Q is the collision integral. The macroscopic parameters like number density n, bulk velocity u, temperature T, and heat flux vector q of metal vapor can be obtained by integrating the velocity distribution function,
respectively. Where m is the mass of metal atom. The expressions of collision integral Q in the BGK model and the S-model are written as
where Pr is the Prandlt number, ω is the collision frequency, FM is local Maxwell distribution,
Non-dimensional quantities are used throughout this article for simplifying the numerical calculation and obtaining common results. Some reference values are introduced to normalize the actual quantities as follows:
where, Kns, λS, vS, and vT are the Knudsen number, mean free path, most probable speed, and average thermal speed at the liquid metal pool, respectively. d is the collision diameter of metal atom. The hard sphere model in collision is employed in both kinetic model equations and the DSMC method for comparison. It should be noted that Kns is a description of collision strength that is suitable for researching arbitrary type of metal.
Substitution of the reference values in Eq. (9) and Eq. (10) into Eqs. (1), (6), (7), (8) firstly and then deleting the superscript ‘*’. The kinetic model equations, collision integral and Maxwell distribution read as,
respectively. Macroscopic parameters in Eqs. (2)–(5) may further be shown as follows:
In order to economize computer storage and cost, two reduced velocity distribution functions g and h are introduced,
the kinetic model equations (11)–(14) are firstly multiplied by 1 and respectively and then integrated over vy in the whole velocity space,
Here, X and Y are the reference values. The macroscopic parameters used in the kinetic model equations are obtained by integrating the new velocity distribution functions g and h,
The stationary kinetic model equations (20) and (21) are solved by finite difference methods. Spatial derivations are replaced by one-order upwind difference which is usually used in transport problem. What is more, nonuniform grids in phase space and velocity space are constructed to match the macroscopic flow properties. Denser grids are adopted for fierce variation of flow properties in the area close to the chamber surface, especially the evaporation source. The velocity space mesh is divided according to the gradients of macroscopic velocity and temperature for precise numerical integration which are used to obtain the macroscopic parameters. Finally, the model equations are iterated constantly for a stable flow field and macroscopic parameters.
4. Results and discussion
4.1. Metal evaporation without collimators
Figures 2(a)–2(d) present the dimensionless results of macroscopic parameters which are calculated by the DSMC method, S-model equation, and BGK model equation under the circumstance of Kns = 1.0. One can see that the results based on the three methods have a good agreement quantitatively in most areas of the evaporation chamber. Furthermore, no obvious difference occurs between two kinetic model equations. Metal atoms collide with one another immediately above the evaporation source when they evaporate from the source surface because of large number density. As a result, a fraction of the metal atoms develops downwards velocities that move them back to the liquid metal pool. Internal energy in metal vapor converts into kinetic energy in the meantime. Accordingly, a rapid expansion flow occurs near the evaporation source. The number density and temperature decrease sharply and the bulk velocity, including the x direction and z direction, increase rapidly. Thereafter, the metal vapor flow reaches the free-molecular or near free-molecular regime as the number density decrease at the area far from the evaporation source. Variations of all macroscopic parameters go slowly. Finally, all of the metal atoms are desublimated on the substrate or walls.
Fig. 2. (color online) Distribution of macroscopic parameters without collimator: (a) n, (b) T, (c) ux, and (d) uz. Left side DSMC, right side kinetic model equations: solid line BGK, dashed line S-model.
Figures 3(a)–3(c) show the variations of n/Kns, T, uz at different Kns over z axis at x = 0.0 cm. Comparison between the kinetic model equations and the DSMC method shows good agreement at the whole line from the center of the evaporation source to the substrate for big Kns (Kns ≥ 0.1). However, discrepancies which exceed the numerical uncertainty appear when Kns is small (Kns ˂ 0.1), especially near the metal pool surface. The discrepancies can be explained as follows: the kinetic model equations are the approximate expression of Boltzmann equation, in which collision integral is treated as a “pulling” to the equilibrium velocity distribution simply. The “pulling” strength is collision frequency. However, deviation is brought in the process of simplification because the change rate of velocity distribution function f(v) comes from three parts according to the Boltzmann equation: particle movement, exogenic action and the collision integral. The fraction of change rate from the collision is positively associated with the number density and the deviation from the equilibrium state in the kinetic model equations. Discrepancy becomes more obvious in the domain where number density is high or velocity distribution is consequently in a nonequilibrium state. The velocity distribution function f(v) is just in accordance with these condition when Kns is small near the evaporation source. As a result, discrepancies occur near the source surface and then transmit downstream along the flow line shown in Fig. 3, although the number density is lower and velocity distribution is closer to an equilibrium state. However, the kinetic model equations can provide reliable results with relatively less computational costs when the Knudsen number is big at the source surface.
Fig. 3. (color online) Macroscopic parameters along the z axis without a collimator at x = 0.0 cm: (a) n/Kns, (b) T, (c) uz. Lines DSMC, crosses BGK, circles S-model; green Kns = 0.01, black Kns = 0.10, blue Kns = 1.00, red Kns = 10.0.
4.2. Metal evaporation with collimators
Figures 4(a)–4(d) show the dimensionless macroscopic parameters contours with a collimator under the circumstance of Kns = 1.0. Where the dimensionless temperatures and the absorption rates of walls below the dashed line are all 0.5. It can be seen that the kinetic model equations provide correct results as above while slight discrepancies of number density that exceed the calculation precision occur above the collimator aperture. In this evaporation condition, a fraction of metal atoms is condensed on the collimators, lower walls and metal ingot when they escape from the liquid metal pool and meet the chamber surfaces mentioned before. The rest reflect diffusely back and impede the diffusion toward walls and substrate. Thus, there are a decrease of uz and a destabilization of velocity distribution embodied in the increase of temperature near the collimator aperture. A similar phenomenon also occurs near the lower walls and ingot surface. When the metal vapor passes through the collimator aperture, the vapor flow is similar to a cold atomic beam because of the low number density and temperature. Finally, a plating spot with clear edge appears at the center of the substrate shown in Fig. 6. It should be noted that the temperature and bulk velocity counted in the DSMC model in the domain above the collimators are not correct because of low number density.
Fig. 4. (color online) Distribution of macroscopic parameters with collimator: (a) n, (b) T, (c) ux, (d) uz. Left side DSMC, right side kinetic model equations: solid line BGK, dashed line S-model.
Fig. 5. (color online) Macroscopic parameters along the z axis with a collimator at x = 0.0 cm: (a) n/Kns, (b) T, (c) uz. Lines DSMC, crosses BGK, circles S-model; black Kns = 0.10, blue Kns = 1.00, red Kns = 10.0.
The sources of discrepancies mentioned before are analyzed as follows according to the discussion of flow properties. The flow of metal vapor is disturbed by walls below the collimators because of the finite absorption rates. Temperature and bulk velocity of metal vapor, that flies to the boundaries, are not equal to those reflected back, especially at the domain where bulk velocity is perpendicular to the boundary. However, the velocity distribution at the collimator aperture is not influenced. Both two types of velocity distributions of vapor join together and then a new nonequilibrium state occurs below the collimators (By contrast, the velocity distribution in the ideal condition is closer to an equilibrium state due to lack of reflection at the same domain). As a result, the discrepancies that are originated from the nonequilibrium state transmit downstream from the collimators to the substrate, especially for small Kns exhibited by Fig. 6. Moreover, we can see that the kinetic model equations are more sensitive to Kns with collimators by the comparison between Fig. 3 and Fig. 5. The kinetic model equations are applied in an opposite smaller scope inside.
Fig. 6. (color online) The z-direction atomic flux along the collimators. Lines DSMC, crosses BGK, circles S-model; black Kns = 0.10, blue Kns = 1.00, red Kns = 10.0.
Fig. 7. (color online) The z-direction atomic flux along the substrate. Lines DSMC, crosses BGK, circles S-model; black Kns = 0.10, blue Kns = 1.00, red Kns = 10.0.
The z-direction atomic flux, which is the product of number density and z-direction bulk velocity, along the collimators and substrate is presented in Fig. 6 and Fig. 7, respectively. The z-direction atomic flux along substrate presents the deposition thickness which is one of the evaluation criteria of the atomic beam. It can be seen that the kinetic model equations provide inaccurate results for small Kns. Even though for big Kns, z-direction atomic flux along substrate calculated by kinetic model equations is not convincing because the deposition thickness profile over substrate is not performed. Actually, the deviation accumulates gradually as metal vapor diffuses from liquid metal pool to substrate. As a result, the deposition thickness profiles over substrate obtained from kinetic methods are not exact even for lower evaporation rate.
5. Conclusion
The metal evaporation of different structures was investigated by the kinetic model equations and the DSMC method. The dimensionless distributions of macroscopic parameters at a series of evaporation rates obtained from three methods were compared. It was shown that whether the collimators were in consideration or not, the results obtained from the kinetic model equations were in good agreement with those from the DSMC method when the evaporation rate was low correspondingly. For the case with a high evaporation rate, the kinetic model equations could only provide correct results qualitatively. And the inaccuracy increased further if the collimators were included in calculations for the diffusion effects. In addition, we found that the profiles of deposition thickness over the substrate from the kinetic model equations were not correct for the beam case. It is believed that the results set up a basis for research into more practical structures and more complex evaporation conditions.
